## Replicate the analysis in Where’s the Pivot? Obstruction and Lawmaking in the Pre-cloture Senate,
## American Journal of Political Science, Vol. 48, No. 4, October 2004, Pp. 758–774

library(MASS)
library(dplyr)

## Figure 1

obstruct <- matrix(scan("obstruction.incidents.txt"),byrow=T,ncol=3)

pdf(file="obstruct.pdf",height=5,width=6.5)

matplot(obstruct[,1],obstruct[,2:3],xlab="Congress",ylab="Number of pieces of legislation",lty=c(1,2),type="l",lab=c(10,20,4),xaxt="n")	

par(cex=.75)
	
cong <- c(1,5,10,15,20,25,30,35,40,45,50,55,60,64)

axis(1, at=cong,labels=c("1","5","10","15","20","25","30","35","40","45","50","55","60","64"))

leg=c("Significantly obstructed","Significantly obstructed and failed")

legend(1,18,leg,lty=c(1,2))

dev.off()

## Table 1

load("pivot_table1_2_ajps04.Rdata")

## pcty: percent voting yea (999999 indicates missing value)
## lt67: =1 if % voting yea < 67
## lt60: =1 if % voting yea < 60
## senate: senate roll call (as opposed to House)
## petersen: legislation indicated as significant by Petersen
## precloture: =1 if vote occurred prior to adoption of cloture rule
## endcong: vote took place within 30 days of end of session
## endcongpre: vote took place within 30 days of end of session in pre-cloture congress
## switsnpr: presidential regime switch, pre-cloture era
## switsnpo: presidential regime switch, post-cloture era
## majpct2: majority party seat share
## lopsided: lopsided vote, defined as > 2/3rds in support

## Results for Table 1

subset(pivot_table1_2_ajps04, pcty != 999999 & senate == 1 & petersen == 1) %>% group_by(precloture) %>% summarise(mean(pcty),mean(lt67), mean(lt60))

## Robustness checks mentioned in text

## average coalition sizes before and at end of congress 

subset(pivot_table1_2_ajps04, pcty != 999999 & senate == 1 & petersen == 1) %>% group_by(precloture,endcong) %>% summarise(mean(pcty),mean(lt67), mean(lt60))

## keeping voice votes and treating them as unanimous or near unanimous

subset(pivot_table1_2_ajps04, pctyvv != 999999 & senate == 1 & petersen == 1) %>% group_by(precloture) %>% summarise(mean(pctyvv),mean(lt67yvv), mean(lopsided))

## Table 2 

table2.df <- subset(pivot_table1_2_ajps04, pcty != 999999 & senate == 1 & petersen == 1)

table2.m1 <- lm(pcty ~ precloture + switsnpr + switsnpo, data=table2.df)

summary(table2.m1)

table2.m2 <- lm(pcty ~ precloture + switsnpr + switsnpo + majpct2, data=table2.df)

summary(table2.m2)

table2.m3 <- lm(pcty ~ precloture + switsnpr + switsnpo + endcongpre, data=table2.df)

summary(table2.m3)

## Robustness checks mentioned in text

## logits with lopsided vote as outcome

lopsided.out1 <- glm(lopsided ~ precloture + switsnpr + switsnpo, data=table2.df)

summary(lopsided.out1)

lopsided.out2 <- glm(lopsided ~ precloture + switsnpr + switsnpo + majpct2, data=table2.df)

summary(lopsided.out2)

lopsided.out3 <- glm(lopsided ~ precloture + switsnpr + switsnpo + endcongpre, data=table2.df)

summary(lopsided.out3)

## Results using either sweep 1 or sweep 2 from Petersen data; NB:
## fn. 27 of the paper implies that we included in the sample
## legislation identified as significant in both sweep 1 and sweep 2,
## when in reality we included in the sample legislation identified as
## significant in either sweep.

sweep.df <- subset(pivot_table1_2_ajps04, pcty != 999999 & senate == 1 & petersen == 1 & (sweep == 1 | sweep == 2))

sweep.m1 <- lm(pcty ~ precloture + switsnpr + switsnpo, data=sweep.df)

summary(sweep.m1)

sweep.m2 <- lm(pcty ~ precloture + switsnpr + switsnpo + majpct2, data=sweep.df)

summary(sweep.m2)

sweep.m3 <- lm(pcty ~ precloture + switsnpr + switsnpo + endcongpre, data=sweep.df)

summary(sweep.m3)

## Results using legislation identified as significant by either
## Petersen or Dell and Stathis; NB: fn. 27 implies that we included
## in the sample for a robustness check legislation identified as
## significant by both Petersen and by Dell and Stathis; instead, the
## robustness check includeds legislation identified as significant by
## either Petersen or Dell and Stathis

ps.df <- subset(pivot_table1_2_ajps04, pcty != 999999 & senate == 1 & (petersen == 1 | stathis == 1))

ps.m1 <- lm(pcty ~ precloture + switsnpr + switsnpo, data=ps.df)

summary(ps.m1)

ps.m2 <- lm(pcty ~ precloture + switsnpr + switsnpo + majpct2, data=ps.df)

summary(ps.m2)

ps.m3 <- lm(pcty ~ precloture + switsnpr + switsnpo + endcongpre, data=ps.df)

summary(ps.m3)

## Table 3
           
load("pivot_table3_ajps04.Rdata")

## run probit with passage as dependent variable
## covariates are as follows:
## coalsize: Coalition size
## daystilendcong: Days left in congress
## coaldaystilend: Coalition size × Days left in congress
## endcongfil: End of congress
## coaldaystilendcong: Coalition size × Days left in congress × End of congress

glm.table3.out <- glm(pass~coalsize + daystilendcong + coaldaystilend + endcongfil + coaldaystilendcong,data=table3.df,binomial(link="probit"))

summary(glm.table3.out)

xint <- cbind(1,table3.df$coalsize, table3.df$daystilendcong, table3.df$coaldaystilend, table3.df$endcongfil, table3.df$coaldaystilendcong)

phat <- pnorm(xint%*%coef(glm.table3.out))

## compute expected percent correctly predicted

epcp <- (1/nrow(xint))*(sum(table3.df$pass*phat + (1-table3.df$pass)*(1-phat)))

cat("epcp: ",epcp,"\n")

## produce predicted probabilities with confidence intervals; NB: the
## values produced here will vary some from what is reported in the
## article because of the use of a random number generator. We did not
## set a seed for the random number generator in the original analysis
## and thus we cannot reproduce those original values. But the values
## produced with this code should be qualitatively the same.

r <- 5000

b.tilde <- mvrnorm(r,coef(glm.table3.out),vcov(glm.table3.out))

xvars <- select(table3.df,c(coalsize, daystilendcong, coaldaystilend, endcongfil, coaldaystilendcong))

xvars.meds <- apply(xvars,2, median)

phat.xvars.meds <- pnorm(b.tilde%*%c(1,xvars.meds))

phat.xvars.meds.quants <- quantile(phat.xvars.meds, c(.025, .5, .975)) 

phat.xvars.meds.quants

## predicted probs at end of congress

xvars.meds.endc <- apply(xvars[xvars$endcongfil==1,],2,median)

phat.xvars.meds1 <- pnorm(b.tilde%*%c(1,xvars.meds.endc))

phat.xvars.meds1.quants <- quantile(phat.xvars.meds1, c(.025, .5, .975)) 

phat.xvars.meds1.quants

## last day

xvars.meds.endc["endcongfil"] <- 1

xvars.meds.endc["coaldaystilend"] <- .51

xvars.meds.endc["coaldaystilendcong"] <- .51

phat.xvars.meds.last <- pnorm(b.tilde%*%c(1,xvars.meds.endc))

phat.xvars.meds.last.quants <- quantile(phat.xvars.meds.last, c(.025, .5, .975)) 

phat.xvars.meds.last.quants 

## biggest coalition in sample 

## not end of congress

xvars.max <- apply(xvars,2, max)

xvars.meds.max <- xvars.meds

xvars.meds.endc <- apply(xvars[xvars$endcongfil==1,],2,median)

xvars.meds.max["coalsize"] <- xvars.max["coalsize"] 

xvars.meds.max["coaldaystilend"] <- xvars.meds.max["daystilendcong"]*xvars.max["coalsize"] 

phat.xvars.max <- pnorm(b.tilde%*%c(1,xvars.meds.max))

phat.xvars.max.quants <- quantile(phat.xvars.max, c(.025, .5, .975)) 

phat.xvars.max.quants

## end of congress

xvars.meds.max["endcongfil"] <- 1

xvars.meds.max["coaldaystilend"] <- xvars.meds.endc["coaldaystilend"]*xvars.max["coalsize"] 

xvars.meds.max["coaldaystilendcong"] <- xvars.meds.endc["coaldaystilend"]*xvars.max["coalsize"] 

phat.xvars.max <- pnorm(b.tilde%*%c(1,xvars.meds.max))

phat.xvars.max.quants <- quantile(phat.xvars.max, c(.025, .5, .975)) 

phat.xvars.max.quants
